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Abstract 

We perform high-statistics Monte Carlo simulations of three-dimensional Ising spin-glass models 
on cubic lattices of size L: the ± J (Edwards- Anderson) Ising model for two values of the disorder 
parameter p, p = 0.5 and p = 0.7 (up to L = 28 and L = 20, respectively), and the bond-diluted 
bimodal model for bond-occupation probability ph = 0.45 (up to L = 16). The finite-size behavior 
of the quartic cumulants at the critical point allows us to check very accurately that these models 
belong to the same universality class. Moreover, it allows us to estimate the scaling-correction 
exponent w related to the leading irrelevant operator: uj = 1.0(1). Shorter Monte Carlo simulations 
of the bond-diluted bimodal models at pb = 0.7 and pi, = 0.35 (up to L = 10) and of the Ising spin- 
glass model with Gaussian bond distribution (up to L = 8) also support the existence of a unique 
Ising spin-glass universality class. A careful finite-size analysis of the Monte Carlo data which takes 
into account the analytic and the nonanalytic corrections to scaling allows us to obtain precise and 
reliable estimates of the critical exponents v and ??: we obtain z/ = 2.45(15) and r/ = —0.375(10). 

PACS numbers: 75.10.Nr, 64.60.Fr, 75.40. Cx, 05.10.Ln 
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I. INTRODUCTION 



The Ising model with random ferromagnetic and antiferromagnetic couphngs is a sim- 
phfied model^ for disordered uniaxial magnetic materials which show glassy behavior in 
some region of their phase diagram, such as Fei_a;Mna;Ti03 and Eui-a^Baa^MnOs; see, e.g., 
Refs. 2,3,4. The random nature of the short-ranged interactions is mimicked by nearest- 
neighbor random bonds. This model is also interesting per se, since it provides a laboratory 
in which the combined effect of quenched disorder and frustration can be studied. 

It is now well established that three-dimensional Ising spin-glass models present a high- 
temperature paramagnetic phase and, for some values of the parameters, a low-temperature 
glassy phase (if the frustration is small, the low-temperature phase is ferromagnetic). The 
two phases arc separated by a continuous phase transition, which is expected to have 
universal features, i.e. to belong to a universality class which is independent of the de- 
tails of the model and, in particular, of the disorder distribution. Several numerical 

^Qj.]^g5,6,7,8,9,10,ll,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33 J^^^g addreSSed theSe iSSUeS, 

considering various Ising spin-glass models, characterized by different disorder distribu- 
tions, with or without dilution. Over the years many estimates of the critical exponents 
have been obtained. We mention the most recent ones for the correlation-length expo- 
nent u: ly ^ 2.39(5),=^° u = 2.72(8),29 ^ ^ 1.5{3),^^ v = 1.35(10),22 v = 2.15(15),2i 

V — 1.8(2),^° obtained from simulations of the symmetric model with bimodal distribution; 

V — 2.22(15) for the bond-diluted symmetric bimodal model with = 0.45;^^ v = 2.44(9)^° 
and u = 2.00(15)^^ for the symmetric model with Gaussian disorder distribution; u = 2.4(6) 
for the random-anisotropy Heisenberg model in the limit of large anisotropy,^^ which is ex- 
pected to be in the same Ising spin-glass universahty class. ^^'^^'^^ Moreover, the analysis of 
different quantities has often given different estimates of the same critical exponent, even 
in the same model. For instance, recent Monte Carlo (MC) studies^^'^'^ find significant dis- 
crepancies among the estimates of the exponent u obtained from the finite-size scaling (FSS) 
at Tc of the temperature derivative of ^/L, of the Binder cumulant, and of the overlap sus- 
ceptibility. For the bimodal Ising model Ref. 30 quotes u — 2.39(5), u — 2.79(11), and 
u = 1.527(8), from the analysis of these three quantities, respectively. A likely reason for 
these discrepancies is the presence of scaling corrections, which may be quite important in 
spin-glass systems since the absence of fast MC algorithms makes it necessary to work with 
systems of relatively small size. 

One of the aims of the present paper is a detailed analysis of the role of scaling corrections 
in spin-glass systems. We present a general renormalization-group (RG) analysis based on 
the Wegner expansion, which allows us to predict the corrections to the asymptotic critical 
behavior for the different quantities. In particular, we show that the analytic dependence 
of the relevant scaling fields on the model parameters, such as the temperature, may give 
rise to scaling corrections that decay as powers of where L is the linear size of the 

lattice. Since u ^ 2.45 in Ising spin-glass systems, they decay quite slowly and may give rise 
to systematic deviations which are difficult to detect, given the small interval of values of L 
which can be considered in MC simulations. Their presence explains some inconsistencies in 
the standard analyses of MC data reported in the literature. Thus, it is crucial to take scaling 
corrections into account for an accurate study of the critical behavior, for a robust check of 
universality among different models, and for reliable estimates of universal quantities such 
as the critical exponents. 

In this paper we report a high-statistics MC study of different Ising spin-glass models. 
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We consider the ±J Ising model for two values of the disorder parameter, the bond-diluted 
symmetric bimodal model with various values of the dilution, and also the model with 
Gaussian disorder distribution. We determine the FSS behavior of several RG invariant 
quantities, such as the ratio = ^/L is the second-moment correlation length) and 
the quartic cumulants defined from the overlap variables. We verify with good precision 
their independence on the model and disorder distribution, providing an accurate evidence 
of universality. Then, we obtain an estimate of the leading correction-to-scaling exponent uu: 
00 = 1.0(1). Finally, we determine the critical exponents. We analyze the MC data at the 
critical point and in the high-temperature phase, taking into account the RG predictions for 
the scaling corrections and the precise above- reported estimate of cu. We obtain 



In this work we extend the results presented in Ref. 32. First, we have significantly increased 
the statistics of the large-L data for the bimodal symmetric Ising model and we have added 
some data for other diluted models and for the bimodal model with Gaussian distributed 
couplings. Second, we present a much more detailed analysis of the critical-point data and 
a new analysis of the high-temperature data obtained in the parallel-tempering simulations. 
This allows us to check the universality of the critical behavior of the correlation length 
and of the susceptibility in the high-temperature phase. Finally, we discuss the extended- 
scaling scheme,^^'^^ which is inspired by the high-temperature expansion. As already noted 
in Ref. 29, this scheme shows an apparent improvement of the scaling behavior with respect 
to the naive approach in which scaling corrections are simply neglected, at least for some 
quantities, e.g., the overlap susceptibility. However, as we shall show, such an improvement 
is only marginal for the purpose of obtaining accurate results. Indeed, this requires to take 
into account the analytic and nonanalytic scaling corrections predicted by the RG theory. 

The paper is organized as follows. In Sec. II we define the models we investigate and 
the quantities that are computed in the MC simulation. In Sec. Ill we derive the FSS 
predictions of the RG theory which are the basis of our FSS analyses. Some details are 
reported in App. A. In Sec. IV and in App. B we give some details on the MC simulations. 
In Sec. VA we discuss universality, verifying that the infinite- volume limit of the quartic 
cumulants and of = ^/L is independent of the model. In Sec. VB we compute the leading 
correction-to-scaling exponent u. In Sec. VI we compute the critical exponents u and r] and 
the critical temperature for the different models by using the data close to the critical point. 
In Sec. VII we present a global analysis of all available high-temperature data obtained in 
our parallel-tempering MC simulations. We again determine the critical exponents and show 
that the FSS behavior of several quantities is universal. Moreover, we discuss the extended- 
scaling scheme of Ref. 29. In Sec. VIII we compute the high-temperature zero- momentum 
quartic couplings. Finally, in Sec. IX we present our conclusions. 



i/ = 2.45(15), 
rj = -0.375(10). 



(1) 



Then, using scaling and hyperscaling relations we obtain 



/5 = z/(l + r?)/2 = 0.77(5), 
7 = (2-r/)z/ = 5.8(4), 
a = 2-3i/ = -5.4(5). 



(2) 
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FIG. 1: Phase diagram of the three-dimensional itJ (Edwards-Anderson) Ising model in the T-p 
plane for p > 1/2, i.e., 1 — p < 1/2. The phase diagram is symmetric mider p ^ ^ — P- 

II. ISING SPIN GLASS SYSTEMS 

A. The ±J Edwards- Anderson Ising model and its phase diagram 

We consider the ±J Edwards- Anderson Ising model on a simple cubic lattice of linear 
size L with periodic boundary conditions. The corresponding Hamiltonian is^ 

H = —''^J^ya^ay, (3) 

where ax = ±1, the sum is over the nearest-neighbor lattice sites, and the exchange interac- 
tions Jxy are uncorrelated quenched random variables with probability distribution 

P{Jxy) = p5{Jxy - 1) + (1 - p)5{Jxy + 1). (4) 

The usual bimodal Ising spin glass model, for which [Jxy] = (brackets indicate the average 
over the disorder distribution), corresponds to p = 1/2. For p 7^ 1/2 we have 

[Jxy] = 2p - 1 ^ 0, (5) 

and ferromagnetic (or antiferromagnetic) configurations are energetically favored. Note that 
the free energy and also the correlations of the overlap variables that we shall define below 
are invariant under p —>■ 1 — p and thus we can always assume p > 1/2. 

The T-p phase diagram of the three-dimensional ±J Ising model is sketched in Fig. 1. 
The high-temperature phase is paramagnetic for any p. The nature of the low-temperature 
phase depends on the value of p: it is ferromagnetic for small values of 1 — p, while it is 
glassy with vanishing magnetization for sufficiently large values of 1 — p. The three phases 
are separated by transition lines, which meet at a magnetic-glassy multicritical point (MGP), 
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usually called Nisliimori point, which is located along the so-called Nishimori line' 
defined by the relation {p > 1/2) 



,38,39,40,41,42 



Tn{p) 



2 



(6) 




On the Nishimori line the magnetic and the overlap two-point correlation functions are equal. 

The paramagnetic-ferromagnetic (PF) transition line starts at the Ising transition at 
p — 1 and extends up to the MGP at p = p*. For p — 1 there is the standard Ising transition 
at^^ Tjs — 4.5115248(14). Disorder is a relevant perturbation of the pure three-dimensional 
Ising fixed point. As a consequence, the Ising point p = 1 is a multicritical point^^ with 
crossover exponent = ais, where^^ ais = 0.1096(5) is the Ising specific-heat exponent. The 
critical behavior for any 1 > p > p* belongs to the randomly-dilute Ising (RDI) universality 
class, characterized by the correlation-length critical exponent^^'^^ u — 0.683(2) and by 
the magnetic exponent 77 = 0.036(1). The coordinates of the MGP in the T-p plane are''^ 
T* = 1.6692(3) and p* = 0.76820(4). The multicritical behavior is characterized by^^ the 
thermal exponent u = 1.64(5), the crossover exponent = 1.67(10), and the magnetic (and 
also overlap) exponent r] = —0.114(3). 

The paramagnetic-glassy (PG) transition line starts at the MGP and extends up to p = 
1/2 (actually up to p = 1 — p* = 0.23180(4) due to the symmetry p ^ 1 — p of the phase 
diagram). A reasonable hypothesis is that the PG critical behavior is independent of p, as 
found in mean-field models. Hence, for any p* > p > 1 — p* the PG transition is expected 
to belong to the same universahty class (named ISG in Fig. 1) as that of the bimodal Ising 
spin-glass model at p = 1/2. The critical behavior along this transition line is the main 
subject of this paper. As we shall see, the universality hypothesis is fully confirmed by our 
FSS analyses at p = 0.5 and p = 0.7. 

The nature of the ferromagnetic-glassy (FG) transition line is not clear yet. At fixed p 
the following inequahty holds:^^'^^ 



where the subscripts indicate the temperature of the thermal average, and Ti^{p) is the 
temperature along the Nishimori line, defined in Eq. (6). This relation shows that ferromag- 
netism can only exist in the region p > p* and that the system is maximally magnetized 
along the Nishimori fine. Ref. 50 (see also Refs. 41,51) also argues that the FG transition 
line coincides with the line p = p*, from T = T* to T = 0. Recent numerical investiga- 
tions^^'^^'^^ of the two-dimensional ±J model have shown that this conjecture is not exact, 
though quantitative deviations are small: at T = the critical value Pc where ferromag- 
netism disappears is definitely larger than p*, indicating a reentrant transition line. In three 
dimensions Ref. 55 quotes Pc — 0.778(5), which is slightly larger than p* = 0.76820(4), with 
an associated critical exponent z/ = 1.1(3): it is therefore likely that the conjecture does not 
hold in three dimensions as well. We also mention that a mixed low-temperature phase,^^ 
in which ferromagnetism and glassy order coexist, is found in mean-field models^^ such as 
the infinite-range Sherrington-Kirkpatrick model.^^ Its presence has been confirmed in the 
±J Ising model on a Bethe lattice. However, there is no evidence of a mixed phase in the 
±J Ising model on a cubic lattice^^ and in related models. In particular, the numerical 
results of Ref. 55 show that the onset of the glassy behavior at T = occurs close to the 
point where the ferromagnetic phase disappears, and are consistent with a single transition 



\[{'^x<yy)T]\ < [\{<yx(ry)TMip)\], 



(7) 
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within numerical precision. Nevertheless, the existence of such a mixed phase is still an open 
problem, as discussed in Ref. 58. 



B. Other Ising spin-glass models 

We also consider the bond-diluted bimodal Ising model (BDBIM) defined by Hamiltonian 
(3) with bond probability distribution 

+ {1-P,)5{j,y) . (8) 

A PG transition occurs for sufficiently small values of f — p^, i.e. for > psG- While 
investigations at T = indicate that psG should be identified with the bond-percolation 
point^°'^^ (Pperc — 0.2488126(5) on a simple cubic lattice^^), recent investigations of the 
critical behavior close to the percolation point suggest that psc is larger than Pperc-^^ 
The model can be extended by considering the distribution 

P{Jxy) = Pb [pSiJccy - J) + (1 - p)SiJ,y + J)] + (1 - Pb)S{J,y). (9) 

In this case, for pb > psG {psG may depend on p) we expect a T-p phase diagram similar 
to the one sketched in Fig. 1 for the ± J Ising model without dilution, with a PF and a PG 
transition line meeting at a multicritical point. 

A PG transition is also observed in the random-bond Ising spin-glass model with Gaussian 
bond distribution: 

PiJxy) = (10) 

This transition is expected to be in the same universality class as that of the bimodal Ising 
model.^° 



xyl 



Pb 



xy 



C. Overlap thermodynamic quantities 

In this work we focus on the critical behavior of the overlap parameter 

^ (11) 

where the spins cr^ ^ belong to two independent replicas with the same disorder realization 
{Jxy}- The corresponding correlation function is 

G{x) = [{qoq.)] = (12) 

where the angular and square brackets indicate the thermal average and the quenched average 
over {Jxy}, respectively. We define the susceptibility x a-^id the second-moment correlation 
length ^ as 



X 



J2G{x), (13) 



G{0)- G{p) 
4sin2(p^i,/2) G{p) ' 
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where p = (pmin, 0, 0), pmin = Stt/L, and G{q) is the Fourier transform of G{x). 
We also define the RG invariant quantities 

= i/L, (15) 

a..M. (16) 

where 

i^k^{C£i^f). (18) 

X 

We call them phenomenological couplings and denote them by it! in the following. 

In the high-temperature paramagnetic phase, we also consider the zero-momentum quartic 
couplings 

- (19) 
G22 ^ (20) 

where 

X4 = (K] - 3[/.^]) , (21) 
X22 ^ m - ■ (22) 
The critical limit T — > T+ of the zero-momentum quartic couplings G4 and G22 is universal. 



III. FINITE-SIZE SCALING 

In this section we summarize some basic results concerning FSS, which allow us to un- 
derstand the role of the analytic and nonanalytic scaling corrections. We consider two Ising 
spin-glass systems coupled by an interaction 

hY.^^-hY,a^^af\ (23) 

X X 

where /i is a constant external field. The model is defined on a cubic lattice of linear size L 
with periodic boundary conditions. 

By applying standard RG arguments we expect the disorder-averaged free-energy density 
to be the sum of a regular part and a singular part: 

^(/3, K L) = ^reg(/3, K L) + .Fsing(/3, K L), (24) 

where /? = 1/T. The regular part is expected to depend on L only through exponentially 
small terms, while the singular part encodes the critical behavior. The starting point of FSS 
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is the scaling beliavior of tlie singular part of the free-energy density (see, e.g., Refs. 36,47, 
64,65): 

^,i,g(/5,/i,L) = L-^F{uHLy\utLy\{v,Ly^}), (25) 

where d is the space dimension, Uh and Ut are the scaling fields associated with h and 
the reduced temperature t ~ 1 — /5//5c (their RG dimensions are i/h = {d + 2 — rj)/2 and 
Ut = l/i', respectively), and Vi are irrelevant scaling fields with < 0. At the critical point 
we have Ut{t — 0,h — 0) — and Uh{t — 0,h — 0) — 0, while, generically, we expect 
Vi{t = 0,h = 0) ^ 0. Since j/j < 0, for large L the free energy can be expanded in powers of 
{vjL^*}. Therefore, we can write 

^sing(/3, h, L) = L-^fiuhLy^utLy^) + v^L'"-'^ UuhLy\utLy') + . . . (26) 

where the leading nonanalytic correction-to-scaling exponent cu is related to the RG dimen- 
sion of the leading irrelevant scaling field v^, = vi, cu — —yoj- The scaling fields are 
analytic functions of the system parameters — in particular, of h and t — and are expected 
not to depend on L. Note also that the size L is expected to be an exact scaling field for 
periodic boundary conditions. For a general discussion of these issues, see Ref. 65, Sec. Ill 
of Ref. 66, and references therein. In general, Ut and Uh can be expanded as 

Uh = huh{t) + 0{h'^), Uhit) = ah + ait + 0{t^), (27) 
Ut = ctt + CQ2t^ + C2o/i^ + C2ihH + 0{f, h\ hH), (28) 

where we used the fact that the free energy is symmetric under h ^ —h. In the expansion of 
Uh,t around the critical point h,t = 0, the terms beyond the leading ones give rise to analytic 
scaling corrections. 

The scaling behavior of zero-momentum thermodynamic quantities can be obtained by 
performing appropriate derivatives of J- with respect to h. For instance, for the overlap 
susceptibility at h — we obtain 

= L^-'^Uh{tyg{utLy^) [1 + v^L-^g^{utLy^) + ...]+ g,,^{P). (29) 

h=0 

The function gregiP) represents the contribution of the regular part T^egiP, h,L) of the free- 
energy density and is L independent (apart from exponentially small terms). It gives rise to 
a correction proportional to L''"^. Analogous formulae hold for the 2n-point susceptibilities. 
The FSS of the phenomenological couplings is given by 

R{(3,L) = r{utLy*)+v^r^{utLy^)L-^ + ... 

^ R* + r'{0)cttLy' + ... + c^L-^ + ..., (30) 

where R* = r(0), — v^r^{0), and the second line holds only very close to the critical 
point, for tL^* ^ 1. A proof of Eq. (30) for the phenomenological couplings U4, C/22, and R^ 
is presented in App. A. 

The thermal RG exponent yt = is usually computed from the FSS of the derivative 
R' of a phenomenological coupling it! with respect to /3 at Pc- Using Eq. (30) one obtains 

i?' = — = Ly*{df,ut) [r'iutLy*) + v^L-^rl{utLy*) + •••]. (31) 



dh? 
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One may also consider the derivative x' = dx/dp of tlie susceptibility x- From Eq. (29) we 
obtain 

+2L'-^Uhd0Uhg{utLy^) + • • • + g[,^{(5). (32) 

Note that the second term in the right hand side gives rise to scaling corrections proportional 
to L~y^ = L~^/^ , while the background term g',.eg{P) leads to corrections proportional to 

At T = Tc, setting t = in the above-reported equations, we obtain: 

R^R* + c^L-^ + ..., (33) 

X = cL2-'?(1 + c^L-- + ...), (34) 

i?' = cLi/'^(l + c^L-" + ...), (35) 

X' = cL^-^'+'^'il + c^L-^ + ... + CaL-^'"" + ...)■ (36) 

Note that, unlike the temperature derivative R' of a RG- invariant quantity, x! ^Iso presents 
an L"^/^ scaling correction, due to the analytic dependence on t of the scaling field Uh (for 
this reason we call it analytic correction). Since, as we shall see, in the Ising spin-glass case 
1/u K, 0.4 and uj ~ 1.0, the scaling corrections in %' decay significantly more slowly than 
those occurring in R' . This makes the ratio 

^ ~ L^/"" (37) 
X 

unsuitable for a precise determination of v and explains the significant discrepancies observed 
in Ref. 30. 

Instead of computing the various quantities at fixed Hamiltonian parameters, one can 
also consider FSS keeping a phenomenological couphng R fixed at a given value i?/.^^'^^ This 
means that, for each L, one determines [3f{L, Rf), such that R{P = Pf{L, Rf), L) = Rf, and 
then considers any quantity at = (3f{L, Rf). The value Rf can be specified at will, as long 
as Rf is taken between the high- and low-temperature fixed-point values of R. For Rf ^ R*, 
where R* is defined in Eq. (33), /?/ converges to Pc as 

- /3e ~ (38) 

while ior Rf — R* we have 

f3f-P,^ L-y---. (39) 
Indeed, if Utj{L, Rf) is the value of Ut for /3 = PfiL, Rf), we obtain from Eq. (30) 

UtjLy^ = B{Rf) + vMRf)L-'' + ■■■ (40) 

where, using Eq. (30), 

riB{x)) = X , (41) 
_ rUBix)) 

'^^^^ r'{B{x)) • ^^^^ 
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Now, if Rf = R*, we have B{Rf) = 0, which imphes Utj ~ L^^*^'^, hence Eq. (39). Other- 
wise, B{Rf) is different from zero and we obtain the behavior (38). 

If we now substitute relation (40) into Eqs. (29), (30), (31), and (32), wc obtain the 
expansion of the different quantities at fixed Rf, which we denote by adding a bar: given 
0{13,L), we define 0{L,Rf) = 0[Pf{L, Rf), L]. For Rf = R*, since Utj ~ L'^'-^ we 
reobtain Eqs. (33), (34), (35), and (36), with different coefficients, of course. If Rf ^ R* , 
we must be more careful. The behavior of another phenomenological coupling it!^ does not 
change qualitatively and we still have 

R^{L,Rf)^Rl^c^L-^ + (43) 

where R*^ = ra[B{Rf)] is universal. It depends on Rf and satisfies R^ = -R* for Rf = R*. 
Also x' behaves as it does at fixed T = Tc, i.e. it follows Eq. (36). On the other hand, x ^^'^ 
R! present additional analytic corrections. Indeed, since [see Eq. (27)] 

Uh,f^ah + -B{Rf)L-y' + --- (44) 

(nonanalytic 0{L~'^) corrections have been neglected), Eq. (29) gives 

g{B{Rf))[l + 0{L-)]. (45) 

If Rf ^ R* , B{Rf) 7^ 0, and thus analytic corrections occur. Note that, if Rf is close to R*, 
since B{R*) — 0, we have 

B{Rf) r^Rf- R*. (46) 

Hence, in this case the analytic corrections are small, of order [Rf — R*)L''^/^ . In general, 
corrections of order L~^/^ have amplitudes proportional to {Rf — R*Y . 



Xa{L,Rf) = L 



2-rj 



al + 2^B{Rf)L- 



yt 



+ 0(L-2^*) 



IV. MONTE CARLO SIMULATIONS 

In the MC simulations we employed the Metropolis algorithm, the random-exchange 
method (often called parallel-tempering or multiple Markov-chain method), and multispin 
coding. See App. B for details on their implementation. 

We simulated the ±J Ising model at p = 0.5 for L = 3-14,16,20,24,28, at p = 0.7 
for L = 3-12,14,16,20, and the BDBIM at pb = 0.45 for L = 4-12,14,16. We averaged 
over a large number Ng of disorder samples: Ng ~ 6.4-10^ up to L — 12, Ng/lO^ ~ 
2400, 2800, 1500, 245, 150, 18, respectively for L = 13, 14, 16, 20, 24, 28 in the case of the ± J 
Ising model at p = 0.5. Similar statistics were collected at p = 0.7 (except for L = 20 where 
the statistics were approximately 1/3 of those for p = 0.5), while for the BDBIM statistics 
were smaller (typically, by a factor of two for the small lattices and by a factor of 6 for the 
largest ones). We also considered the BDBIM at pb — 0.7 and pb — 0.35 and the Ising model 
with Gaussian distributed couplings, but only performed simulations for small values of L: 
1 = 4,6, 8, 10 for the BDBIM and L = 4, 5, 6, 8 for the Gaussian model. 

For each L and model we performed parallel-tempering runs. This allowed us to estimate 
the different quantities in a large interval [/?min, /3max]- To fix /?max we used the results of 
Ref. 30, which provided the best estimates of Rt at the time we started our simulations: 
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TABLE I: Quartic cumulant U4 at fixed i?^ = 0.63 for the ± J Ising model at p = 0.5 and p = 0.7, 
for the BDBIM at pb = 0.45, 0.7, 0.35, and for the Ising spin-glass model with Gaussian bond 
distribution. 



L 


±'^p=0.5 


±'^p=0.7 


BDBIMp,=o.45 


BDBIMp,=o.7 


BDBIMp,=o.35 


Gaussian 


4 


1.48231(6) 


1.46813(5) 


1.49036(8) 


1.48480(6) 


1.49164(9) 


1.49145(13) 


5 


1.48985(6) 


1.47597(6) 


1.49853(8) 






1.4996(2) 


6 


1.49446(6) 


1.48193(6) 


1.50300(9) 


1.49618(9) 


1.50788(9) 


1.5033(2) 


7 


1.49753(6) 


1.48642(6) 


1.50544(9) 








8 


1.49987(6) 


1.48984(6) 


1.50714(9) 


1.50082(9) 


1.51320(13) 


1.5063(5) 


9 


1.50136(6) 


1.49260(6) 


1.50815(9) 








10 


1.50273(6) 


1.49478(6) 


1.50889(9) 


1.50382(11) 


1.5146(2) 




1 1 

J. L 


X.OUOOOl L) J 












12 


1.50469(6) 


1.49781(7) 


1.50984(13) 








13 


1.50541(11) 












14 


1.50618(10) 


1.50030(12) 


1.5103(3) 








16 


1.50702(13) 


1.50220(13) 


1.5113(3) 








20 


1.5081(3) 


1.5048(5) 










24 


1.5089(4) 












28 


1.5108(13) 













0.627(4) and 0.635(9) for an Ising model with bimodal and Gaussian distributed bonds, 
respectively. Thus, in most of the runs /3max was chosen so that R^{Praax,L) ~ 0.63. Only 
in the most recent simulations (two runs with L — 20 and 24) of the ±J Ising model at 
p = 0.5 did we take /9max such that R^{Prna.x, L) ~ 0.66. We checked thermalization by using 
the recipe outlined in Ref. 30, see App. B for some details. 

As already emphasized in Refs. 46,69, in high-precision MC studies of random systems 
one should be careful when computing disorder averages of products of thermal expectations, 
for instance the quartic cumulant U22 defined in Eq. (17). Indeed, naive estimators have a 
bias that may become significantly larger than the statistical error if A^<j is large. We use 
(essentially) bias-free estimators, defined as discussed in Ref. 46; some details are given in 
App. B. 

In total, the MC simulations took approximately 40 years of CPU time on a single core 
of a 2.4 GHz AMD Opteron processor. 

V. UNIVERSALITY AND CORRECTION-TO-SCALING EXPONENT oj 

A. Analysis of the quartic cumulants at fixed 

In this section we investigate the universality of the critical behavior of Ising spin-glass 
models by comparing the limit for L ^ 00 of f/4 and U22 computed at fixed i?^, denoted by 
f/4 and U22, respectively. As discussed in Sec. Ill, for sufficiently large L, tJ^ and U22 are 
expected to behave as 

= + c#L-, (47) 
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TABLE II: Quartic cumulant tJ22 at fixed = 0.63 for the ± J Ising model at p = 0.5 and p = 0.7, 
for the BDBIM at p}, = 0.45, 0.7, 0.35, and for the Ising spin-glass model with Gaussian bond 
distribution. 



L 


±Jp=o.5 


±Jp=o.7 


BDBIMp,=o.45 


BDBIMp,=o.7 


BDBIMp,=o.35 


Gaussian 


4 


0.13714(6) 


0.13955(6) 


0.13581(9) 


0.13602(7) 


0.14635(10) 


0.13522(15) 


5 


0.14088(7) 


0.14087(6) 


0.13935(10) 






0.14014(24) 


6 


0.14277(7) 


0.14181(7) 


0.14166(10) 


0.14193(10) 


0.14501(10) 


0.14227(24) 


7 


0.14392(6) 


0.14262(7) 


0.14308(10) 








8 


0.14478(7) 


0.14318(7) 


0.14415(10) 


0.14414(10) 


0.14597(15) 


0.1434(5) 


9 


0.14522(7) 


0.14380(7) 


0.14470(10) 








10 


0.14561(7) 


0.14418(7) 


0.14518(10) 


0.14536(12) 


0.14586(24) 




1 1 

J. L 


vj. -L^oyoi i J 












12 


0.14618(7) 


0.14465(7) 


0.14605(14) 








13 


0.14650(11) 












14 


0.14671(10) 


0.14531(13) 


0.1462(3) 








16 


0.14675(14) 


0.14553(14) 


0.1469(4) 








20 


0.1469(4) 


0.1458(5) 










24 


0.1477(5) 












28 


0.1496(14) 













where the constants are universal (model independent) but depend on the fixed value of 
R^; cu is the leading scaling-correction exponent. 

In Tables I and II we report the estimates of U4 and U22 at fixed — 0.63 for different 
models. Without performing any analysis, one can immediately observe that the results 
obtained for the different models are very close, and appear to approach the same large-L 
limit as L increases. For instance, the estimates of U4 for the largest lattices differ at most 
by 0.5%, while those of U22 vary by a few percent. This already provides strong support to 
universality. 

For a more detailed analysis, let us first consider the three models for which we have most 
data, the ± J model at p = 0.5 and p = 0.7, and the BDBIM at pb = 0.45. The MC estimates 
of Ui{L) are shown in Fig. 2 versus 1/L. The results for the ±J Ising model at p = 0.5 and 
p — 0.7 fall quite nicely on two straight lines that approach the same point as L — > 00. They 
support the universality of the critical behavior and show the presence of scaling corrections 
with an exponent ^ 1.0. In the case of the BDBIM with ph = 0.45, the data apparently 
show a faster approach to the same infinite-volume limit. A fit of (/^{L) to + cL~^ using 
all data with L > 4 gives an effective exponent e ^ 2. However, for sufficiently large L, 
L > 10 say, scaling corrections are consistent with a linear 1/L behavior, as in the other 
models. The corresponding amplitude \ci\ is quite small, at least a factor of two smaller than 
in the undiluted bimodal model. This means that for L < 10 next-to-leading corrections to 
scaling dominate and determine the apparent approach of the data to the infinite-volume 
limit. Since the ratios of the amplitudes of the leading scaling corrections are universal, a 
small IC4I implies that the leading nonanalytic scaling correction is small for any quantity. 
Thus, in this model the approach of any thermodynamic quantity to the critical limit should 
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FIG. 2: (Color online) Estimates of U/i{L) vs L ^, for the it J model ai p = 0.5 and p = 0.7, and 
the BDBIM at ph = 0.45. The dotted lines are drawn to guide the eye. 
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FIG. 3: (Color online) Estimates of U22{L) vs L ^, for the ±J model at p = 0.5 and p = 0.7, and 
the BDBIM at p^ = 0.45. The dotted lines are drawn to guide the eye. 

be faster than in the ±J models without dilution, as already noted in Ref. 28, at least for 
sufficiently large lattice sizes, where next-to-leading corrections can be neglected. Similar 
results are obtained for f/22, see Fig. 3. However, in this case next-to- leading corrections 
appear to be more important, since the 1/L behavior is observed for somewhat larger values 
of L. 

In order to estimate the universal infinite- volume values f/4 and U22, we perform fits of 
f7# to + c^L~'^ . The exponent e is either a free parameter, or is set equal to 1, which 
corresponds, as we discuss below, to our best estimate of the leading scaling-correction 
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TABLE III: Results of the fits of U4 and U22 at fixed = 0.63 to U* + cL~^ for several values 
of Linin, the minimum lattice size allowed in the fits. DOF is the number of degrees of freedom of 
the fit. The fits labelled "p = 0.5 h p = 0.7" are fits in which the data for p = 0.5 and p = 0.7 are 
analyzed together, assuming the universality of [/| and C/I2. 



Model 




-^min 


ut 


xVdof 




xVdof 


± J,p = 0.5 


free e 


6 


1.5127(4) 


0.7 


0.1478(3) 


0.7 






8 


1.5135(9) 


0.8 


0.1487(9) 


0.7 






10 


1.5119(11) 


0.4 


0.1482(12) 


0.8 




£=l 


8 


1.5144(2) 


0.8 


0.1490(2) 


0.6 






12 


1.5139(4) 


0.6 


0.1488(4) 


0.9 


±J,p = 0.7 


free e 


6 


1.5153(7) 


2.0 


0.1481(8) 


0.8 






8 


1.5145(14) 


2.7 


0.1471(9) 


0.9 




e = l 


8 


1.5143(2) 


2.3 


0.1478(2) 


0.9 






12 


1.5153(5) 


0.1 


0.1482(5) 


0.8 


BDBIM,p6 = 0.45 


free e 


6 


1.5120(3) 


0.5 


0.1479(6) 


0.6 






9 


1.5132(22) 


0.6 


0.1485(23) 


0.5 




e = l 


8 


1.5153(3) 


0.8 


0.1496(3) 


0.4 






12 


1.5148(12) 


1.2 


0.1489(13) 


1.0 


± J,p = 0.5 & p = 0.7 


e = 1 


8 


1.5143(1) 


1.2 


0.1485(1) 


2.1 






12 


1.5145(3) 


0.9 


0.1486(3) 


0.8 






16 


1.5133(9) 


0.4 


0.1489(10) 


1.1 




£=1.2 


8 


1.5120(1) 


4.9 


0.1479(1) 


2.5 






12 


1.5127(3) 


0.4 


0.1482(3) 


0.9 






16 


1.5123(8) 


0.3 


0.1486(9) 


1.1 



exponent u. In Table III we report the results. They are all consistent with the estimates 

tJl = 1.514(2), [7*2 = 0.148(1). (48) 

The accuracy and stability of the fits can also be appreciated in Fig. 4, where the fit results 
are plotted versus the minimum lattice size Lmin allowed in the fits. We can thus conclude 
that the estimates of the quartic cumulants for the ±J Ising model and the BDBIM at 
Ph — 0.45 are fully consistent with universality. The relative differences are approximately 
one per mille in the case of and one per cent for U22- 

We also considered the BDBIM for the other dilution values, i.e., for ph = 0.7 and pb = 
0.35, though in this case we have data only up to L = 10. The estimates of are shown 
versus 1/L in Fig. 5. For p^ = 0.7 scahng corrections clearly behave as 1/L and are larger 
than at pb — 0.45. For the scaling-correction amplitude C4 defined in Eq. (47) we obtain 
the estimate C4 ~ —0.10, to be compared with the result C4 ~ —0.11 for the ±J model at 
p = 0.5 and C4 ^ —0.05 for the diluted model at pb = 0.45. For pb = 0.35 we do not have 
enough data to estimate C4; however, the data reported in Fig. 5 apparently approach the 
infinite-volume hmit faster than at Pb — 0.45. As indicated by the MC data for pb < 0.27 of 
Ref. 63, the scaling corrections should increase as Pb further decreases. This suggests that 
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FIG. 4: Estimates of and C/22' obtained by fitting the quartic cumulants C/^ at fixed = 0.63 
to U^ + cL~^. Here L min is the minimum lattice size allowed in the fits. The fits labelled "p — 0.5 
& p = 0.7" are fits in which the data for p = 0.5 and p = 0.7 are analyzed together. The dotted 
lines correspond to the final estimates = 1.514(2) and = 0.148(1). 



the optimal model — the one with the smallest leading scaling corrections — corresponds to a 
dilution in the range 0.3 ^ pi ^ 0.4: the model with pb = 0.35 should be close to the optimal 
one. 

The estimates of U4 and U22 for the Ising spin-glass model with Gaussian bond distribution 
are shown in Figs. 5 and 6: they are very close, and clearly approach, the estimates (48). 
The results for L = 8 differ by 0.5% {U4) and 3% (f/22) from the asymptotic value. 

In conclusion, the results for the ±J model, for the BDBIM, and for the model with 
Gaussian distributed coupUngs provide a very strong evidence of universality for the critical 
behavior of Ising spin-glass models along the PG transition line. 
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FIG. 5: Estimates of U4:{L) for the BDBIM for several values of and for the Ising spin-glass model 
with Gaussian bond distribution. The dotted lines correspond to the estimate (48), i/| = 1.514(2). 
The dashed straight lines associated with the data at pi, = 0.7 and ph = 0.45 show the linear 
behavior of the data for their largest lattices. 
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FIG. 6: [/4 (below) and U22 (above) for the Ising spin-glass model with Gaussian bond distribution. 
The dotted lines correspond to the estimates (48), C7| = 1.514(2) and U22 = 0.148(1). 

B. The RG exponent of the leading irrelevant operator 

The analyses of O4 and U22 also give estimates of uj. The most precise ones are obtained 
from fits of U4. In Fig. 7 we show the estimates of u obtained from fits of 1/4^ to + CpL~'^ 
and from fits of the difference U4{p = 0.5; L) — lf4{p = 0.7; L) to hL~^ . We estimate 

cu= 1.0(1). (49) 

Consistent but less precise results are obtained from fits of U22- The result (49) is consistent 
with the less precise estimates uj = 0.84:~^q^y and cu = 1.0(4) reported in Refs. 21 and 27, 
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FIG. 7: (Color online) Estimates of the leading correction-to-scaling exponent uj versus the min- 
imum lattice size Lmin used in the fit. The fits labelled "p = 0.5 & p = 0.7" are fits in which 
the data for p = 0.5 and p = 0.7 are analyzed together. The dotted lines correspond to the final 
estimate to = 1.0(1). 



respectively. 

In order to verify that scaling corrections are controlled by a single correction-to-scaling 
term with exponent uj = 1 and we have not determined an effective exponent arising from 
several almost degenerate exponents, we check that the ratio C22/C4 is universal, where c# 
is the scaling-correction amplitude appearing in Eq. (47), as predicted by standard RG 
arguments. In order to determine this ratio, we fit together the available estimates of U (L) 
for the ±J Ising models at p = 0.5 and p = 0.7 to U{L) = U* + CpL~'^, fixing 6 = 1 and 



taking U* and the two amplitudes Cp= 
L > L^in = 10 give 



:0.5 



and c„=o.7 as free parameters. Fits of the data for 



f/4* = 1.5142(2), C4(p = 0.5) = -0.114(2), C4(p = 0.7) = -0.194(2), (50) 
(x^/DOF ^ 1.3, where DOF is the number of degrees of freedom of the fit), and 

f/*2 = 0.1484(2), C22(p = 0.5) = -0.027(2), C22(p = 0.7) = -0.044(2), (51) 
(xVDOF ^ 1.4). It follows 

C22/C4 = 0.24(2) for p = 0.5, (52) 

C22/C4 = 0.23(1) for p = 0.7, (53) 

which are in good agreement. These results are quite stable with increasing Lmin- For 
Lmin = 12 we obtain C22/C4 = 0.24(3) and C22/C4 = 0.24(2) for p = 0.5 and p = 0.7, while for 
-^^min = 14 we find C22/C4 = 0.19(9) and C22/C4 = 0.21(5), respectively; in both cases the fits 
have x^/DOF < 1. Moreover, the variation of the ratio C22/C4 with respect to a change of 
the exponent e by ±0.1 [it corresponds to the uncertainty of u, u = 1.0(1)], is smaller than 
the above- reported errors. These results fully support our interpretation of the 0{L~^'^) 
corrections as the corrections arising from the leading irrelevant scaling field. 
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FIG. 8: (Color online) Crossing point Pctoss{L) for the it J Ising model at p = 0.5, from several pairs 
L, 2L, as obtained from and C/4, versus ^-^/'^-'^ ^ L~^'^. The dotted lines show the result of a 
fit to Pc + cL~^'^, which takes into account only the data satisfying L > Lmm = 8 — it corresponds 
to « 0.054. The bar along the y-axis corresponds to the final estimate Pc = 0.902(8). 




FIG. 9: (Color online) Estimates of /3c as obtained from fits of and C/4 to Eq. (56) with iv = 1, 
for the ibJ Ising model at p = 0.5. All fits have x^/DOF < 1.5. The lines correspond to the final 
estimate Pc = 0.902(8). 
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FIG. 10: (Color online) The values of the critical temperature Tc versus 1 — p for the various 
transitions of the itJ Ising model. The estimates at the paramagnetic-ferromagnetic transitions 
are taken from Ref. 43 for the Ising transition, from Ref. 44 for the RDI transition, and from Ref. 42 
for that along the line. 

VI. CRITICAL TEMPERATURE AND EXPONENTS 



A. Estimates of [5c 

In order to determine the critical temperature, we perform a standard FSS analysis of 
and the quartic cumulants. Figure 8 shows the crossing points [3ctoss{,L) determined by 
solving the equation 

R{PcToss, L) = -R(/3cross, 2L), (54) 

for the ± J Ising model a.t p = 0.5, using and f/4. In the large-L limit /5cross(-^^) is expected 
to converge to (3c as 

(3cross{L) - ~ L-i/'^-- ~ (55) 

since u ~ 2.45 and 00 ~ 1.0. We perform a combined fit of R^ and U4 to jSc + c^L'"^ with 
e = —1.4, taking Pc and the two amplitudes cr^ and cu^ as free parameters. Using only 
the resuhs with L > L^in = 8, we obtain (3c = 0.9035(22) with xV^OF ^ 1.1 (since the 
estimates of i?^ and U4, are correlated the error is only indicative). The corresponding lines 
are drawn in Fig. 8. We also consider larger values of Lmin- We find that the estimates of (3c 
vary significantly (much more than the statistical error) when Lmin varies, indicating that 
the systematic error due to the additional scaling corrections is significantly larger than the 
statistical one. Equivalently, one can perform fits to 



R{L,(3c 



R* + bL- 



(56) 



taking R 
of the results on L 
estimates 



R^ or U/^J^ The results are shown in Fig. 9. Taking into account the dependence 



and their variation as cu varies by one error bar we obtain the final 
/?c = 0.902(8), Te = 1.109(10). (57) 
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FIG. 11: (Color online) Estimates of the exponent z/ from the analysis of R'^, at = 0.63 
and R^ = 0.654, for the itJ Ising model at p = 0.5, obtained by fitting the data without scaling 
corrections (wsc), Eq. (59), and with scaling corrections, Eq. (60). We only show the results of the 
fits which satisfy x^/DOF < 1.5, up to values of -Lmin where the errors blow up. The dotted lines 
correspond to the final estimate = 2.45(15). 



We regard the error as conservative. We also consider the pseudo critical /3f{L) defined in 
Sec. Ill at fixed R^, which converges to /3c as L — oo, cf. Eqs. (38) and (39). The results are 
consistent with the estimate (57). The analysis of the crossing points and the fits to Eq. (56) 
also provide estimates of the limit L ^ oo of and f/4 at the critical point. We obtain 

R* = 0.645(15), f/; = 1.50(2). (58) 

Concerning the other models, similar analyses of R^ and U4 give f3c = 0.87(1) [Tc = 
1.149(13)] for the ±J model at p = 0.7, and Pc = 1.54(2) [T^ = 0.649(8)] for the BDBIM at 
Pb = 0.45. The corresponding values of -R| and are consistent with the estimates (58). 

Finally, in Fig. 10 we collect all results for the critical temperature of the ±J Ising 
model at its various PF and PG transitions as a function of the disorder parameter p. 



B. The critical exponent v 

We estimate the critical exponent u from the finite-size behavior of R'^ = dR^/d(3 and 
U'^ = dU^/dp at a fixed value R^j of R^. As we discussed in Sec. Ill, the best choice 
for R^j is -R|, otherwise R'^ and U!^ present analytic corrections. In the present case, the 
estimate (58) of is not very precise, hence the corrections of order L^^^'^ ~ cannot 
be completely neglected. However, since their amplitude is proportional to {R^ — R^), we can 
assume that they are small for R^j close to Thus, in order to estimate the systematic 
error induced by them, we proceed as follows. We repeat the analysis of R' at fixed R^ 
using two different values of R^ and neglecting in both cases the L~^^'^ corrections: we use 
i?^ = 0.630 and R^ = 0.654, which are close to our best estimate = 0.645(15). For 
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FIG. 12: (Color online) Estimates of the exponent u from the analysis of 0^, at = 0.63 and 
= 0.654, for the it J Ising model at p = 0.5, obtained by fits without scaling corrections (wsc), 
Eq. (59), and fits with scaling corrections, Eq. (60). We only show the results of the fits which 
satisfy x^/DOF < 1.5, up to values of L^i^ where the errors blow up. The dotted lines correspond 
to = 2.45(15). 



both values of we determine an estimate of u. Our final result is obtained by linear 
interpolation to = 0.645. The systematic error is estimated from the difference of the 
results at R^ = 0.630 and i?^ = 0.645. 

The exponent u is obtained from fits of R'^ and U!^ to 

Ini?' = a + -InL, (59) 

Ini?' = a + hnL + bL'', (60) 

with e = uj = 1.0(1). The results of the fits of R'^ for the ±J Ising model with p = 0.5 are 
shown in Fig. 11 as a function of the minimum lattice size Lmin allowed in the fits. They are 
quite stable with respect to Lmin, depend very little on the fixed value i?^, and change only 
slightly as u varies by ±0.1, corresponding to one error bar. In particular, the fit of R'^ at 
R^ = 0.654 to Eq. (60) gives z/ = 2.52(4) [2] (z/ = 2.53(7) [1]) for L > L^i^ = 8 (L,^in = 10); 
here the error in brackets takes into account the uncertainty on uj. Analogously, the fit of 
the data at i?^ = 0.630 gives z/ = 2.51(4) [2] (z/ = 2.47(7) [1]) for the same values of Lmm- In 
both cases xVDOF ^ 1.1 (xVDOF ^ 1.3). Therefore, we obtain z/ = 2.52(6) and 2.51(8) 
at -Rg = 0.645 for Lmin = 8, 10. The systematic error due to the analytic corrections is 
small, 0.01 and 0.04 for the two values of Lmin- The results from fits of f/4, which are shown 
in Fig. 12, favor a smaller value of u, although in substantial agreement. Indeed, fits with 
L^in = 10 give u = 2.30(9) at R^ = 0.654 and z/ = 2.42(9) at R^ = 0.630, thus suggesting 
the estimate z/ = 2.35(9). The error due to the analytic corrections is apparently larger than 
that obtained in the analysis of R'^, approximately 0.07. As our final result we take 

z/ = 2.45(15), (61) 
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FIG. 13: (Color online) Estimates of the exponent v from the analysis of B!^ and 0^, at i?^ = 0.63 
for the ibJ Ising model at p = 0.7, obtained by fitting the data without scaling corrections (wsc) 
and with scaling corrections, Eq. (60). We only show the results of the fits with x^/DOF < 1.5, up 
to values of Lmin where the errors blow up. The dotted lines correspond v = 2.45(15). 



which is consistent with the results obtained for tj'/^ and B!^. 

Results at i?g = 0.63 for the ±J model at p = 0.7 and the BDBIM at pf, = 0.45, shown in 
Figs. 13 and 14, respectively, are fully consistent with the estimate (61). For the ±J model at 
p = 0.7, fits with e = 1.0 and Lmin = 8 of and U'^ at fixed = 0.63 give u = 2.49(5) and 
u = 2.37(6), respectively. In the case of the BDBIM at pb = 0.45 and again for = 0.63, 
fits without scaling corrections are in full agreement. Fits of R'^ to (60) assuming e = 2.0 
are also consistent: for Lmin = 7 they give u = 2.52(4). On the other hand, the fits of 
with e = 2.0 give results that are systematically larger; for instance, we obtain u = 2.61(5) 
for Lmin = 7. This slight discrepancy is probably due to the fact that scaling corrections 
with exponent u = 1 are small in the BDBIM at pb = 0.45 but not completely negligible: 
thus, the fits assume that corrections vanish faster than they really do, obtaining a slightly 
incorrect asymptotic estimate. 



C. The critical exponent rj 

We estimate the exponent r] by analyzing the susceptibility x fixed R^. Its critical 
behavior is discussed in Sec. Ill: at fixed R^, it behaves as 

X{L, R^) = aL^-'' [I + a„(i?5 - i?*)L-i/^ + ■ ■ ■ + a^L"- + . . . + 0^^-^+'' + ..■], (62) 

where the 0(L^^^^) term is the background contribution, cf. Eq. (29). Since R^ is not exactly 
equal to we must take into account the 0{L~^/'') corrections. As for z/, the systematic 
error they induce is estimated by comparing the results at i?^ = 0.630 and i?^ = 0.654, which 
are close to the best estimate Rt = 0.645(15). We perform fits with and without scaling 
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FIG. 14: (Color online) Estimates of the exponent u from the analysis of R'^ and U'^, at i?^ = 0.63 
for the BDBIM at ph = 0.45, obtained by fitting the data without scaling corrections (wsc) and with 
scaling corrections, Eq. (60). We only show the results of the fits with reasonable values of x^/DOF, 
up to values of Lmin where the errors blow up. The dotted lines correspond to = 2.45(15). 
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FIG. 15: (Color online) Estimates of the exponent r/ from the analysis of x{L, R^) for R^ = 0.63 
and R^ = 0.654, for the itJ Ising model at p = 0.5. Only results of fits with x^/DOF < 1.5 are 
shown. The dotted lines correspond to the estimate r/ = —0.375(10). 



corrections, i.e. to 



InX = c + {2 — r7)lnL, 

InX = c + {2 — ■?7)lnL -|- c^U 



(63) 
(64) 
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with e = 1 (the leading scahng correction arising from irrelevant scaling fields) and e = 
2.4 2 — rj (to check the relevance of the background term, which might be large for small 
lattices) at = 0.63 and = 0.654. Fig. 15 shows the estimates of r] for the ±J model at 
p = 0.5. As already mentioned, the systematic error due to the neglected L~^f^ corrections 
is estimated from the difference of the estimates at R^ = 0.63 and R^ = 0.654. Following 
the same method used to determine we obtain r] — — 0.375(2){8}, where the first error 
is the statistical one while the second error takes into account the systematic error due to 
the residual 0{L~^/'^) corrections and corresponds to the uncertainty on our estimate of 
Slightly smaller but perfectly consistent results are obtained in the analyses of the data 
for the other models. For example, in the case of the ±J Ising model at p = 0.7, a fit of 
the data at i?^ = 0.63 to Eq. (64) with e = 1.0 gives r] = -0.381(5) for L > L^in = 10 
(xVDOF ^ 1.6). In the case of the BDBIM at pb = 0.45, a fit of the data at R^ = 0.63 to 
Eq. (64) with e = 2.0 gives rj = -0.385(4) for L > L^in = 9 (xVDOF ^ 1.4). Taking also 
these results into account, we arrive at the final estimate 

rj = -0.375(10). (65) 

We finally discuss the behavior of the derivative x' = dx/ d/3 of the susceptibihty, which in 
some numerical works, see, e.g. Ref. 30, has led to contradictory results for the exponent v. 

We show in the following that this discrepancy is essentially due to the analytic 0{L~^^^) 
corrections discussed in Sec. Ill, cf. Eq. (36). At fixed R^, x' is expected to behave as 

R^) = bL'' [1 + baL-^/" + ■■■ + b^L-'' + ... + bbL"" + •••], (66) 

where a = 1/ v + 2 ~ rj. Using our best estimates of v and r], wc obtain a = 2.78(4). Unlike 
the case of x and the derivative of the phenomenological quantities, the amplitude ba of the 
0{L~^/'^) corrections does not vanish at Tc and thus at R^ — R^; therefore, it is not expected 
to be small when 7?^ ^ The 0{L~^) term in Eq. (66) is the background contribution, 
cf. Eq. (32). In Fig. 16 we show the estimates of a obtained in fits of Inx' to 

Inx' = c + crlnL + c^L"^, 

\nx' ^c + a\nL + ciL-'' +C2L-'\ (67) 

for several values of s, 81,82- The results are consistent with the expected value a ~ 2.78, 
when the 0{L~^/^) corrections are taken into account. 

VII. HIGH-TEMPERATURE ESTIMATES 

In the parallel-tempering simulations we have collected a large amount of data for several 
values of (3 in the high-temperature phase. They have not been used in the analyses we 
have presented in the previous Sections. Here, we shall present analyses that consider all 
high-temperature results. They fully confirm the critical-point estimates discussed above 
and provide additional evidence that all models belong to the same universality class. 
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FIG. 16: (Color online) Estimates of the critical exponent a = \/v+2 — r], as obtained by analyzing 
x' at Rc = 0.654. The dotted lines correspond to fi = 2.78(4). 
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FIG. 17: (Color online) Plot of i/L vs {(3c - 0)1'^/" with (3c = 0.902 and v = 2.45. Data for the 
lb J model with p = 0.5. 



A. Finite-size scaling behavior of ^ and estimates of v 



We determine v from the FSS behavior of the correlation length. The starting point is 
the FSS equation 

^ = f{u,L^'l + "-j^UutL'^n + ■■■, (68) 

where f{x) and fuj{x) are universal (once the normalization of the scaling fields has been 
fixed) and have a regular expansion in powers of x. In order to ensure a finite infinite- volume 
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limit, for x oo they behave as 



f^{x) ~ X 



u{ui-l) 



(69) 



The scahng field Ut is an analytic function of (5 that vanishes at the critical point. Hence, it 
has an expansion 

ut = [^,-f5 + 6(/9c - + 0[(/5e - f5f]. (70) 

In Fig. 17 we report ^{(3, L)/ L versus [jSc — (3)L^^^ using the data for the ±J model at p = 0.5 
and the estimates of j3c and v obtained in the previous Sections: j3c = 0.902 and u = 2.45. 
The data collapse quite precisely onto a single curve. This indicates that the nonanalytic 
scaling corrections are small and so are the analytic ones due to the nontrivial dependence 
of Ut on p. 

In order to obtain quantitative estimates of u, we follow Ref. 71 and perform three different 
fits of ^{(3, L) I L. In the first fit we neglect the nonanalytic scaling corrections, set Ut = Pc~P, 
and fit the data to (fit a) 



Pnix) 



x={i3,-i5)L^I\ 



(71) 



where P„(a;) is a polynomial of degree n. We have chosen polynomials for computational 
convenience, but any other complete set of functions can be used. Note that the specific 
choice (71) of fitting function satisfies the large-x behavior (69). In this fit, the free parame- 
ters are the (n+ 1) coefficients of the polynomial Pn{x), the critical inverse temperature f3c, 
and uJ"^ The critical-point value is equal to Pn{x = 0)^*^/". 

In order to take into account the analytic corrections (fit b), we slightly modify the Ansatz 
(71), setting x = /3c — P + HPc — I^Y and taking h as an additional free parameter. Finally, 
we consider the nonanalj^ic scahng corrections. We use the Ansatz (fit c) 



- —v/n 



Pn(x) + —(l + axr''Qn(x) 



x=(/3,-/3)lV-, (72) 



where Pn{x) and Qn{x) are both polynomials of degree n, and a is a free parameter. We can 
check that Ansatz (72) has the correct infinite- volume limit. For L — >■ oo at fixed j3, i.e., for 

X 



k that Ansatz {72) has the correct mfanite-volu 
oo, we have P„(a;) ~ PnX^, Qn{x) ~ qnX^i and 

r 1 1 -'^Z" 



PnX" + —[axY''qnx'' 



(73) 



where A = ujv. In these fits n must be chosen so that P„ provides an accurate parametriza- 
tion of the scaling function. We have tried several values of n, until the of the fit does 
not change significantly as n increases by 1. In practice, we have taken n between 6 and 10. 
In fit (72) corrections to scaling are parametrized by a polynomial that has the same degree 
as that parametrizing the leading behavior. Our data are not so precise and corrections are 
not so large to require such a large number of parameters. To reduce the number of free 
parameters we have taken n = 6, 9 and set 



Qn{x) = go + ^3^;^ H \- (InX"' 



(74) 
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FIG. 18: (Color online) Estimates of u for the it J model at p = 0.5. We report estimates obtained 
by fitting ^/L to the Ansatz (71) for several values of Lmin and /3min- The dotted lines correspond 
to the estimate u = 2.45(15). 

This choice, which is completely arbitrary, significantly reduces the number of free parame- 
ters, but still allows us to parametrize accurately (at the level of the statistical errors) the 
scaling corrections. 

In order to take into account the additional scaling corrections which are not taken into 
account by our fit Ansatze, we have repeated each fit several times; each time we only include 
data such that /? > /?min and L > L^in- In Fig. 18 we report the estimates of z/ obtained in 
fits of ^/L to Eq. (71) for the ±J model at p = 0.5. Corrections to scaling are clearly visible 
for small values of Pram and L^m- Indeed, at fixed Lmin the estimates decrease, becoming 
approximately independent of /5min for /5min ^ 0.65. The dependence on Lmin is instead tiny, 
and all results with Lmin > 10 are consistent within errors. 

The results of fits that also take into account the analytic and the nonanalytic corrections 
(with a; = 1) are reported in Fig. 19. As far as ly is concerned, no significant differences 
are observed and in all cases the results become independent of Pram for Pmia ^ 0.65. All 
results (with their statistical errors) lie in the interval 2.3 < v < 2.5, and are therefore in 
perfect agreement with the estimate v = 2.45(15) obtained before. Corrections to scaling 
are more evident in the analyses of Pc- Indeed, while analyses without nonanalytic scaling 
corrections give estimates of Pc that cluster around 0.895, those that take the corrections 
into account give values which are somewhat larger. In any case, all results are consistent 
with the estimate Pc = 0.902(8) obtained in Sec. VIA. Finally, these analyses also provide 
estimates of Analyses without nonanalytic scaling corrections give i?| = 0.632(5), while 
those which include scaling corrections give a somewhat higher value i?| = 0.648(5). Also 
in this case, scaling corrections appear to be quite relevant. These results are perfectly 
consistent with the estimate (58), = 0.645(15). The analyses that take into account the 
analytic scaling corrections (fit b) also give estimates of the constant b that appears in the 
expansion (70) of Ut: they vary somewhat with P^ia and give approximately < 6 < 0.3. 

In order to verify universality, we compute u and determine the FSS curves also for the 
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FIG. 19: (Color online) Estimates of v (top) and /3c (bottom) for the it J model aX, p = 0.5. We 
report estimates obtained by fitting ^/L to Eq. (71) (fit a), to Eq. (71) with x = [{(3c — /?) + t>{Pc — 
Pf]L^/'' (fit b), to Eq. (72) with w = 1 (fit c), and to the extended-scahng Ansatz (86) (fit d). The 
dotted lines correspond to the estimates v = 2.45(15) (top) and (3c = 0.902(8) (bottom). 

±J model at p = 0.7 and the BDBIM at pt, = 0.45. We report here only estimates of u 
obtained by fitting the data to Eq. (71), since we do not have data precise enough to allow 
for a detailed study of the scaling corrections. In any case, the results for p = 0.5 indicate 
that scaling corrections do not play much role in the determinations of i^. For the ± J model 
at p = 0.7 we obtain: z/ = 2.382(5) [z/ = 2.427(10)] for /J^in = 0.59 and L^m = 7 [L^m = 9]; 
z/ = 2.347(10) [u = 2.382(20)] for /?mm = 0.68 and the same values of L^in. For the BDBIM 
at pf, = 0.45 we obtain: z/ = 2.574(6) [u = 2.637(7)] for /5min = 0.82 and again Lmin = 7 
[^min = 9]; u = 2.409(12) [u = 2.437(18)] for /^^m = 1-02. These results are in very good 
agreement with the estimate u = 2.45(15) reported above. 
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FIG. 20: (Color online) Plot of ^/L vs {pc - P)L^^'' for the ±J model at p = 0.7 (top) and for the 
BDBIM at pb = 0.45 (bottom). We use v = 2.45, f3c = 0.87 (top), and Pc = 1-54 (bottom). The 
dashed line is the curve R^{cx), where R^^ix) is the function (76), which is estimated in fits of ^/L 
for the lb J model at p = 0.5, and c is a model-dependent constant: c ~ 1.026 for the it J model at 
p = 0.7 and c « 0.5641 for the BDBIM at pb = 0.45. 



In order to verify the universality of the FSS curves we first fitted the data for the ±J 
model at p = 0.5 presented in Fig. 17. Taking u = 2.45, Pc = 0.902, and using only the data 
satisfying L > 10, P > 0.62, we obtain 

| = i?5(x) x = (/3,-/?)Ll/^ (75) 

with (this expression is valid in the interval of values of x for which we have data, i.e., for 
0<x<1.5) 

R^{x) = (6.2828 + 16.8612a; - 39.3317a;2 + 1926.5102a;3 
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FIG. 21: (Color online) Plot of U4 vs ^/L for the ±J model at p = 0.5 and 0.7 and the BDBIM 
at pb = 0.45. 

-17659.3388x^ + 88711.1141a;^ - 256918.2481a;^ + 446776.5137x^ 
-452723.8074x^ + 243001. 4960a;^ - 50040.5243x^°) . (76) 

The function Rs^{x) is universal apart from a rescaUng of its argument. Thus, if we plot ^/L 
vs X = {Pc — P)L^^'^ in any model, the data should fall on the curve R^{cx), where c is a 
model-dependent constant. In Fig. 20 we report the data for the ±J model at p = 0.7 and 
the BDBIM at ph = 0.45. The results show very good scaling and fall on top of the curve 
computed from the data of the bimodal model at p = 0.5, confirming universality. 

We finally consider the cumulant t/4. In the critical limit U4 should be a universal function 
of ^/L, independent of the model. Corrections scale as L~^h{^/L), where h{x) is also 
universal, apart from a multiplicative constant. The numerical estimates of f/4 and ^/L are 
reported in Fig. 21. We consider the ±J model sX p = 0.5 and p = 0.7 and the BDBIM at 
Pd = 0.45. We only consider the largest lattices, so that nonanalytic scaling corrections are 
not visible on the scale of the figure (a detailed study of the L""^ corrections for ^/L = 0.63 
is reported in Sec. VA). All points fall quite precisely onto a single curve, confirming that 
the PG transition in these three models belongs to the same universality class. 



B. Finite-size scaling of the susceptibility and estimates of r] 



We now turn to the determination of the critical exponent 77. The starting point is 
Eq. (29), which we rewrite as 



x{P,L) = L'-^u,{Pfg{utL'/'') 



. Vuj{P) , j\lv\ I 



(77) 



This equation is not very convenient since it involves u^, hence the critical temperature, and 
the exponent v. To reduce the number of unknown parameters, we note that Eq. (68) can 
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FIG. 22: (Color online) Plot of x^'^ (upper panel) and of xi^ 



(lower panel) vs ^/L with 



T] = —0.375. The dashed line in the upper panel is the universal curve C{^/L), as estimated in fits 
of X to Eq. (80). Data for the it J model at p = 0.5. 



be inverted to give 



UtL 



F{^/L) + H^F^{^/L) + 



Inserting in Eq. (77) we obtain the scaling form 



(78) 



(79) 



In Eq. (79) we have singled out and ^""^ instead of L^"'' and L~'^ . With this choice 
C{x) and Cu){x) are regular for x — 0. Note the presence of the function Uh{f3). In the 
FSS limit ^ oo, L oo, at fixed C,/L, we always have (3 — >■ /3c, so that asymptotically it 
should be possible to replace Uh{P) with the constant Uh{Pc)- Therefore, this function gives 
rise to scaling corrections, that we have named analytic corrections in the previous sections. 
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FIG. 23: (Color online) Estimates of rj for the it J model at p = 0.5. We report estimates obtained 
by fitting x to the Ansatz (80) (fit a) and (81) (fit b) for several values of Ly^m and /3min- The 
dotted lines correspond to the estimate rj = —0.375(10). 

In order to understand their relevance for our data, in Fig. 22 (upper panel) we plot C~'^x 
versus ^/L for the ±J model at p = 0.5. It is evident that the data do not fall onto a 
single curve: the scaling-field term Uh{l3) varies significantly with (3 and therefore cannot be 
neglected. 

The previous discussion indicates that, in order to estimate accurately the exponent r/, it 
is essential to include the analytic corrections in the fitting function. We perform two fits. 
In the first one (fit a) we neglect the nonanalytic scaling corrections and consider 

In ^ = -r^ In e + Pn{i/L) + g^(/3), (80) 

where Pn{x) and Qmi.^) are polynomials of degree n and m, respectively. Moreover, we 
require Qm(0) = in order to avoid the presence of two constant terms. As before, n and m 
are varied till the quality of the fit does not change significantly by varying the parameters 
by 1. In practice, we take 6 < n, m < 10. To include the scaling corrections, we also consider 

In ^ = -r^lne + Pn{^lL) + Q^{(3) + C'^S.ii/L), (81) 

where Sp{x) is a polynomial of degree p: we take p < 3. Note that here, as we already did 
in the analysis of ^/L, we neglect the /? dependence of the scaling field v^^. 

The results of the fits for the ±J model at p = 0.5 are reported in Fig. 23. The estimates 
of the fit to Eq. (80) are very stable and show a very tiny dependence on /5min and Lmin. 
For instance, for /3mm = 0.55 we obtain rj = —0.3636(8) (Lmin = 8), r] = —0.3659(27) 
(Lmin = 14), while for /3mm = 0.75 we have r] = -0.3666(10) (Lmm = 8), r] = -0.3657(31) 
(Lmin = 14). Fits with nonanalytic scaling corrections are less stable. We observe significant 
fluctuations that indicate that the data are not precise enough to be sensitive to this type of 
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FIG. 24: (Color online) Scaling-field function UhiP) for the it J model at p = 0.5, as determined 
in the fits of x to the Ansatz (80), for different values of /^min and L^ain- It is normalized such 
that UhiPc) = 1- We also report the approximation UhiP) ~ iP/ Pc)^'^^"^ , which is used in the 
extended-scaling scheme (see the discussion in Sec. VII C). 

scaling corrections. This is consistent with what is observed at the critical point: While the 
results for rj depend strongly on the chosen value for R^j, indicating that the analytic scaling 
corrections are important, essentially no dependence is observed on the nonanalytic ones; 
see, e.g.. Fig. 15. In any case, all results are consistent with the estimate r] = —0.375(10) 
obtained in Sec. VI C. 

The fits also give estimates of the function Uh{l3) that appears in Eqs. (77) and (79). In 
Fig. 24 we plot Uh{,l3) = Uh{/3)/uh{/3c) [uhiP) is normalized so that UhiPc) = 1] as obtained 
in the different fits. The results corresponding to different values of L^in and Pram agree 
nicely, supporting the scaling Ansatz (79). A simple expression which reproduces the results 
reported in Fig. 24 is 

UhiP) = 1 + 0.556247(1 - /3/0.902) + 1.83322(1 - /3/0.902)^ (82) 

which is valid for 0.55 < (3 < 0.902. Once UhiP) has been determined, we can compute the 
scaling function C{^/L) = Uh{l3c)^C{^/L) by considering L''~^x£t/j(/5)~^. Such a combina- 
tion is shown in Fig. 22 (lower panel): all points fall on top of each other, confirming the 
validity of the FSS Ansatz. Moreover, as expected, we find that C(0) is finite and C{^/L) is 
approximately constant for ^/L < 0.15, two properties which are not obvious from the upper 
panel of Fig. 22. These conclusions are consistent with the FSS results for x(2L, /5)/x(L, /3) 
(in this quantity the analytic function Uh{P) cancels out) reported in Refs. 20,28, which 
show that this ratio has a tiny dependence on ^/L up to ^/L ^ 0.15 (for ^/L = 0.15 we 
have xi'^L, j3)/x{L, f3) ~ 1.02). Note that the curve in the lower panel of Fig. 22 (which 
corresponds to the dashed line in the upper panel) is the limiting curve of the points that 
appear in the upper panel as L — >■ oo. Since the rate of convergence is very slow (at fixed ^/L 
data converge as L~^^'^), it is clear that such an asymptotic behavior can only be observed 
on enormously large lattices. Thus, in order to estimate r], it is crucial to take the function 
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Uh{P) into account. 

The function C{x) is universal, apart from a model-dependent multiplicative constant. 
We write it as 

C{x) = br{x) r(0) = 1, (83) 

where r{x) is universal. A fit of the data reported in Fig. 22 gives 

r(a;) = 1 + 5.9622y - 104.4625y^ + 1516.2443y^ 

-12260.6638y^ + 50105.6104y^ - 80471.2150/, 

y = exp{-l/x), (84) 

and b ^ 1.9395. The expression (84) is valid for x = ^/L < ^ 0.645. 

The same analyses can be repeated for the ±J model ai p = 0.7 and the BDBIM at 
Pb — 0.45. Here we only present results corresponding to fits to the Ansatz (80). Our data 
are not precise enough to allow us to perform fits which include the nonanalytic scaling 
corrections. The results are consistent with the estimate 1] = —0.375(10). For the ±J model 
at p = 0.7, we obtain r] = -0.366(2) {p^i^ = 0.59) and r] = -0.366(3) (/5min = 0.68) for 
Lmin = 7, and r] ^ -0.368(2) {(3min = 0.59) and 7] = -0.366(3) {/J^in = 0.68) for L^in = 9. 
For the BDBIM, we obtain 77 = -0.361(2) and -0.359(3) for L^in = 7 and 9 and any P^nin 
in the range [0.82, 1.12]. 

As a further check of universality we determine the scaling behavior of the combination 
Xi^~'^^h^- III Fig- 25 we report this quantity for the ±J model at p = 0.7 (upper panel) 
and the BDBIM at = 0.45 (middle panel), using in both cases r] — —0.375 and our best 
estimate of Uh- As expected all points fall onto a single curve. If universality holds, these 
curves should be parametrized as bT{x)^ where r(a;) is given in Eq. (84) and 6 is a model- 
dependent constant. In the case of the ± J model (upper panel of Fig. 25) we observe a very 
good agreement, while for the BDBIM at pb — 0.45 (middle panel) some discrepancies occur 
for ^/L > 0.4. There are two reasons for them. First, the function Uh is not precisely known 
for (3 ~ I3c- in this range of values of (3 it varies slightly (10%) with /9min and Lmin- Second, 
the plot depends on r]. If we use 77 = —0.360, which is the value obtained in the fits for the 
diluted model which provide Uh, we obtain the lower panel of Fig. 25. Discrepancies are now 
significantly reduced. 

It is worth noting that the functions Uh are approximately the same in the three models 
we study, if one considers them as a function of the reduced temperature t = 1 — j3 / j3c- 
For < t < 0.4 — this is the interval of t which is probed by our simulations — the ratio 
{t) I Uh,raoA<i\2{t) is coustaut, withiu our precision, for any pair of models. This result 
is somewhat unexpected within RG theory, because these functions are not universal. 

Finally, let us comment on the FSS approach of Ref. 73, applied to spin-glass systems in 
Refs. 20,28. In this approach one considers the ratio x(2L, /3)/x(L, /3). This choice has a 
significant advantage: the scaling-field function Uh{l3) cancels out, so that the leading scaling 
corrections are the nonanalytic ones. As we have shown here, they are quite small, so that 
very good scaling is observed and reliable infinite-volume estimates are obtained. Analytic 
scaling corrections come in again when considering the critical limit of the infinite- volume 
results Xoo(/5)- Indeed, since A = uju ^ 2.45, for jS ^ I3c the analytic corrections dominate: 

Xoo(/3) = (/?c - /?)-^ [&0 + &i(/9c - /?) + &2(/3c - (3f + feA(/?c -(5)^ + ■■■]. (85) 
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FIG. 25: (Color online) Plot of xC'^'^u^^ vs i/L for the ±J model at p = 0.7 with = -0.375 
(upper panel) and for the BDBIM at pi, = 0.45 with r/ = —0.375 (middle panel) and r/ = —0.360 
(lower panel). In each panel we report the curve bT{^/L), where T{^/L) is defined in Eq. (84) and 
6 is a model-dependent constant. We use b = 1.99346, 2.12794, 2.13923, in the upper, middle, and 
lower panel, respectively. They are determined by requiring a perfect fit for ^/L < 0.2. 
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FIG. 26: (Color online) Extended-scaling results for the it J model at p = 0.5. In the upper panel 
we report vs ^/L. The dashed line is the universal curve C{^/L), as estimated in fits of 

X to Eq. (80). In the lower panel we report estimates of ry obtained in three different fits for several 
values of Lmin and /3mm- Fit c uses the Ansatz (89), fit d and e are defined in text, below Eq. (89). 
The lines correspond to the estimate rj = —0.375(10). 

C. Extended-scaling scheme 

In this Section we consider the extended-scahng scheme introduced in Ref. 29. It consists 
in a particular choice of scahng variables, which, according to Ref. 29, should somehow 
decrease scaling corrections and thus allow a faster convergence to the critical limit. Let us 
consider first C,/L. In this scheme the appropriate fit Ansatz is 

= P„(x)-'^/", X ^ {Pi - /5^)(L//?)l/^ (86) 
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where Pn{x) is a polynomial of degree n. The results for the ±J model at p = 0.5 are 
reported in Fig. 19 and should be compared with those obtained by fits to Eq. (71), which 
neglect any scaling corrections (fit a). They are substantially equivalent and no improvement 
is observed. This can be explained by noting that, for the present values of /?c and i/, we 
have 

(/?2 - /?2)/3-V- ^ i.88(/?, - /3)[1 - 0.10(/?e -/?) + ...]. (87) 

Thus, fit (86) is essentially equivalent to a fit with analytic corrections (fit b) with b = —0.10. 
Such a value of b is small — hence, this change of the scaling variable does not have much 
influence on the final results — and is close to what wc obtain numerically, though not fully 
consistent (we predict < 6 < 0.3). We also tried fits with x = (/?^ — /3^)L^/'', which 
should be the natural variable in spin-glass systems, given the symmetry under /? —>■ — /?.^^ 
These fits are significantly worse than the previous ones for /3 < 0.70. For larger values, no 
significant differences arc observed. These results can be understood by noting (/3^ — /?^) = 
1.80(/3c — f3)[l — 0.55(/3c — (3) + . . .]. Thus, this choice of scaling variable corresponds to 
assuming b = —0.55 in Eq. (70), which is significantly larger than what we find numerically. 
Therefore, if we use {(3^ — (3^) as approximate thermal scaling field, the analytic corrections — 
in this case they are proportional to (/9^— /3^)^ — are more important than in the case in which 
Ut is simply approximated with Pc — P- 

The extended-scaling scheme can also be applied to the analysis of the susceptibility. It 
amounts to consider the scaling Ansatz^^ 

x{f3,L)=f3^-'e-'C{C/L). (88) 

In Fig. 26 wc show versus ^/L for the ±J model at p = 0.5. Scaling is better 

than that observed in the upper panel of Fig 22: the scatter of the data points is signifi- 
cantly reduced, indicating that (3'^^^ approximates the scaling-field term uf^ better than a 
constant. However, the rescaled data are still far from the asymptotic curve (the dashed 
line) determined numerically above, indicating that the residual analytic scaling corrections 
are also in this case not negligible. This is better understood, by comparing the function Uh, 
as determined in the fits, with the approximation {(3 / (3c)^^'^~^ , which follows from Eq. (88). 
As can be seen from Fig. 24, the approximate expression proposed in Ref. 29 has the correct 
qualitative shape, but differs significantly from the quantitative point of view. For these 
reasons, wc do not expect the scaling Ansatz (88) to be particularly useful to estimate r] 
from our data. 

To understand the role of the residual analytic scaling corrections on the determinations 
of 77, we fit the data for the ± J model at p = 0.5 to the scaling Ansatz (fit c) 

ln^ = -r7ln| + P„(e/L), (89) 

where Pn{x) is a polynomial of degree n. The results are reported in Fig. 26. They vary 
strongly with j3„ii^ and Lmin, indicating that scaling corrections are sizable and not negligible. 
As a test we have also repeated the fits to Eqs. (80) and (81) replacing Inx/^^ with lnx/9^/^^ 
and In^ with ln^//3 (we call fit d and fit e the fits corresponding to Eqs. (80) and (81), 
respectively). As expected, the results are identical to those obtained in fit b and fit c, 
respectively. Indeed, the fits only differ in the parametrization of the analytic function 
Uh- Note that the extended-scaling approximation for Uh is not analytic in f3, since f3 = 
is a branching point. This is, however, irrelevant in practice, since we are looking for 
approximations of Uh in the interval 0.6 < /3 < 0.9, which is quite far from (3 — 0. 
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VIII. THE ZERO-MOMENTUM QUARTIC COUPLINGS G4 AND G22 IN THE 
HIGH-TEMPERATURE PHASE 

In this section we consider the zero-momentum quartic couphngs G4 and G22 defined in 
Eq. (19) and (20) and estimate their infinite-volume critical value defined by 

hm hm G#{L,P). (90) 

We consider the ±J Ising model at p = 0.5 and use the estimates of the quartic couplings for 
(3 = 0.55 and L = 12, 14, 16, 18, and 32, and for (3 = 0.625 and L = 16, 18, 28, 32, and 40. 
We combine results obtained in the random-exchange runs that we performed for our FSS 
study around /3c and results obtained in standard MC simulations for these two values of 
pJ"^ First, we investigate the infinitc-vohimc limit. The correlation length converges rapidly 
(see Fig. 27). For instance, for P = 0.55 we obtain ^ = 1.7888(3), 1.7885(5), 1.7872(3) for 
L = 18, 20, 32, while for /? = 0.625 we have ^ = 3.2694(12), 3.2709(13) for L = 32 and L = 
40: for L/^ > 10 the results vary by less than 0.1%, indicating that the difference from their 
thermodynamic limit is within 0.1%. Thus, we can take as infinitc-vohimc estimates those 
obtained on the largest lattices. The quartic couplings show larger finite-size corrections. As 
shown in Figs. 27 the infinite- volume limit is approximately reached for L/^ > 12, within 
our statistical precision. Indeed at /3 = 0.625 we find G4 = 90.46(9), 90.27(16), 90.23(29) 
and G22 = -16.53(5), -16.36(9), -16.40(9) for L = 18, 20, 32, respectively. We can thus 
take the estimates on the largest lattices as infinite-volume estimates. 
In the critical limit we expect 

G#(L^oc,P)^G*^ + c#r^, (91) 

with Lu = 1.0(1). The comparison of the results at /3 = 0.55 and (3 = 0.625 leads us to the 
estimates 

= 90.3(5), G^^^-l^l), (92) 

where the error takes also into account the effects of the 0{^~'^) scahng corrections, which 
is roughly estimated from the difference of the infinite-volume results at the two values of 
f3. These results can be compared with those obtained in Ref. 27 for the random-anisotropy 
Heisenberg model in the limit of infinite anisotropy: G^ = 88(8) and G22 ~ There is 

a substantial agreement, which confirms that the critical behavior of the random-anisotropy 
Heisenberg model for infinite anisotropy belongs to the Ising spin-glass universality class. ^^'^^ 



IX. CONCLUSIONS 



In this paper we discuss the critical behavior of three-dimensional Ising spin-glass systems 
with the purpose of verifying universality, clarifying the role of scaling corrections, and 
determining the critical exponents. More precisely, our results can be summarized as follows. 

i) By using the RG we derive the behavior for L 00 and j3 (3c oi several quantities 
which are routinely measured in MC simulations. In particular, we show that the 
analytic dependence of the scaling fields on the model parameters may give rise to 
corrections which behave If they are neglected, FSS analyses give 
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FIG. 27: (Color online) Estimates of ^(L,/5) (top), Gi{L,^) (middle), and G22{L,l3) (bottom) 
versus L/^ for (3 = 0.55 and /? = 0.625, corresponding to the infinite- volume correlation lengths 
^oo ~ 1-79 and ^oo ~ 3.27. Results for the ±J model at p = 0.5. The dotted lines correspond to 
the estimates = 90.3(5) (middle panel), G22 = —17(1) (lower panel). 
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inconsistent results. Tliese corrections have been overlooked in previous FSS studies. 
Note that the general expressions we obtain are relevant also in other glassy systems, 
in which u is typically large. 

ii) We determine the leading nonanalytic correction-to-scaling exponent uj. We obtain 
ijj — 1.0(1). Note that in Ising spin-glass systems nonanalytic scaling corrections decay 
faster than in the Ising model, in which 00 ^ 0.8, see Ref. 47. The exponent uj is also 
significantly larger than that at the PF transition, which occurs for small frustration: 
UJ = 0.29(2).^^ 

iii) We accurately verify universality. A careful analysis of the quartic cumulants U4, and 
U22 at fixed = 0.63 shows that their limit for L — > 00 is independent of the model 
and of the disorder distribution. The results obtained in the different models differ 
at most by approximately one per mille in the case of U4 and one per cent for 1/22- 
They support the existence of a unique Ising spin-glass universality class. Universality 
is also supported by the FSS analyses of ^ and x iii the high-temperature phase. We 
verify that the FSS curves for these two quantities are independent of the model. 

iv) We determine the critical exponents. For this purpose we perform analyses at the 
critical point and analyses which take into account all high-temperature data. Results 
are consistent, once the analytic and the nonanalytic corrections are taken into account. 
Moreover, they do not depend on the model and disorder distribution. Again, this 
supports the universality of the paramagnetic-glassy transition. We obtain 

i/ = 2.45(15), = -0.375(10). (93) 

Using scaling and hyperscaling relations, we find /3 = i/(l + r))/2 = 0.77(5), 7 = 
(2 - r])iy = 5.8(4), and a = 2 - 3i/ = -5.4(5). 

Our estimates of the critical exponents can be compared with those reported in the literature. 
Earlier estimates of u are reported in the introduction and in Table I of Ref. 30. Some of the 
most recent ones are close to our final estimate. For the exponent r], we quote here the most 
recent results: r] = -0.395(17), -0.37(5),^° r] = -0.40(4),29 and rj = -0.349(18). They 
are all in substantial agreement with our result, which is however significantly more precise. 
We also mention the estimate f3 = 0.52(9), obtained in Ref. 26 by an out- of- equilibrium 
simulation. 

The estimates (93) slightly differ from, and have larger errors than, those obtained in 
Ref. 32: u = 2.53(8) and rj = —0.384(9). There are two reasons for that. First, we have 
significantly extended the runs for L = 20, 24 for the ± J model at p = 0.5. The new results 
have slightly shifted the estimates of /3c, and of the critical exponents. Second, we have 
been more conservative. With our present error bars, the estimates are fully consistent with 
the results of all analyses for the three models we considered. 

We also analyzed our data by using the extended-scaling scheme proposed in Ref. 29. This 
approach might partly take into account the scaling corrections arising from the analytic de- 
pendence of the scaling fields on the reduced temperature but it neglects the nonanalytic 
corrections arising from the irrelevant operators. In some cases, for instance for the over- 
lap susceptibility, this scheme shows an apparent improvement of the scaling behavior with 
respect to the naive approach in which the analytic corrections are simply neglected. How- 
ever, the approximate expressions which follow from the extended-scaling scheme are not 
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sufficiently precise for a liigli-precision study of tiie critical behavior: if one aims at accurate 
estimates, it is necessary to determine the corrections directly from the data. 



APPENDIX A: FINITE-SIZE BEHAVIOR OF THE PHENOMENOLOGICAL 
COUPLINGS 

We now provide a detailed proof of Eq. (30) for the phenomenological couplings U4, U22, 
and = ^/L. First, we present a simple physical argument that clarifies the origin of 
the function Uh{t) introduced in Eq. (27), and then, a more formal argument based on the 
Wegner expansion. 

A short and physical proof goes as follows. First, note that x is not RG invariant and, in 
particular, it is not invariant under field redefinitions. Given a theory with fields 0, define 



X<p = J d'xMx)m)]- (Al) 
If we change variables (p — Zip, we obtain 

X4>-Z' J d'x[{iP{x)m)] = Z\^- (A2) 

Thus, under a field redefinition x ^ Z'^X- The function Uh{t) is exactly the field redefinition 
that relates the bare lattice fields to the renormalized fields. Hence, we expect at criticality 

X = ulL'-^. (A3) 

A t dependent prefactor is absent in the phenomenological couplings, because these quantities 
are RG invariant and in particular, they do not depend on the normalization of the fields. 
This is quite obvious from their definitions. In the continuum limit they can be written as 

TT _ J dxidx2dx3dx4[{(j)ixi)(f){x2)(t>{x3)(f){x4:))] , , 

U4 = 2 ) 1-^4 j 

[/ dxidx2{(f){xi)(l){x2))] 

jj ^ / dxidx2dx3dxi[{(t){xx)(t){x2)) (0(x3)0(a:4))] _ ^ ^^^^ 

[/ dxidx2{(t){xi)(l){x2))\^ 

^2 ^ 1 / dXidX2{Xi - X2Y[{(t){Xi)(l>{x2))] 

6 / dxxdx2{(t){x\)(t){x2)) 

The formula we report for ^ correspond to Eq. (14) for L — > 00, disregarding corrections of 
order L~^. Since the number of fields in the denominator is equal to the number of fields in 
the numerator in all expressions, each quantity is invariant under field redefinitions. 

This discussion clarifies the origin of the t-dcpendent prefactor in Eq. (A3) and is at the 
basis of the statement, which is well accepted in the literature, that Ui{Tc), U22{Tc), and 
^{Tc)/L approach universal constants as L — > cxo. The presence of an analytic prefactor 
would violate universality. Indeed, imagine that a phenomenological coupling R behaves as 

R{(5, L) - K,,(i)^°'^^P°"^''t/(MtL^/'^) + scaling corrections (A7) 
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close to the critical point. Since Uh{t) is model dependent, at the critical point R{pc, L) 
would not be universal for L — > oo. Hence, on purely physical grounds, there must be no 
i-dependent prefactor in any phenomenological coupling, hence no analytic corrections. 

This short discussion shows that and x/L^"'' are conceptually two very different ob- 
jects. The first quantity is RG invariant and shows a universal FSS behavior. The second 
quantity is not RG invariant and, for instance, x{Tc) / L'^~'^ converges to a model- dependent 
constant as L — >■ cxo. 

Eq. (30) can also be derived from the usual Wegner's scaling expression for the free energy. 
Let us first consider U4. Using Eq. (26), we find 



^L^y'^ulf^^\o,utLy*) + --- (A8) 

h=0 

= L''y'^ulf^^\0, UtV') + • • • (A9) 



dh? 



h=0 



[the dots correspond to nonanalytic scaling corrections and bulk contributions, and the 
derivatives refer to the first variable appearing in the scaling function f{x,y)] so that 

which proves Eq. (30). 

To discuss U22 one should generalize Wegner's scaling expression (see Sec. 3.1 of Ref. 46 
for a detailed discussion). Define Z{f3,h,L) as the partition function of two systems at 
inverse temperature (3 defined on a lattice of size coupled by an interaction 

h^cri^cr2x- (All) 

X 

Then, consider 

J^{(3, hi, h2, L) = L-^ [In Z(/3, hi, L) In Z{P, h^, L)] . (A12) 

A scaling Ansatz like Eq. (25) allows one to obtain an expression analogous to that obtained 
for C/4 and to prove Eq. (30) for C/22- 

In order to determine the scaling behavior of we consider a momentum-dependent 
magnetic field. The argument goes as follows: Define Z{f3, h, L,p) as the partition function 
of two systems at inverse temperature defined on a lattice of size coupled by an interac- 
tion h o'ix(J2x cos(p • x). Then, consider the corresponding disorder- averaged free-energy 
density 

J^iP, h, L, p) = L-'' [In Z(/3, h, L,p)]. (A13) 

Under RG transformations L — > XL, momenta scale as p — > p/X, so that the singular part 
of the free-energy density scales as 

T^^^{l3,h,L,p) = L-''f{pL,Uh{h,t,p)Ly\ut{h,t,p)Ly'), (A14) 

where we have neglected the nonanalytic scaling corrections and now the scaling fields depend 
also on p. Taking derivatives with respect to h and then setting /i = 0, we obtain for the 
two-point function [of course Uh{h,t, —p) — Uh{h,t,p)] 

G{p) = u,{t,pfL''-^f^''\pL,0,Ut{0,t,p)Ly^), (A15) 
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where we write as before 

uhih,t,p)^huhit,p) + 0{h^), (A16) 

and we have neglected subleading terms. For p ^ 0, because of the cubic symmetry of the 
lattice, we have 

Uh{t,p)=Uh{t)+0{p^), (A17) 
ut{0,t,p) =ut{0,t) + O{p^), (A18) 

where Uh{t) and Ut{0,t) are the usual (zero-momentum) scaling fields. Hence, for p ^ 0, 
disregarding corrections of order p^, we can express G{p) in terms of the scaling fields that 
appear for p = 0: 

G{p) = UH(t)'L'-^f^^ (pL, 0, ut(0, t)Ly*) + O(p') + • • • (A19) 

In the definition (14) of the correlation length ^ we should consider p = g ~ 1/L. Thus, 
disregarding terms of order we have 

G(0) - G(p) /P'(O.O.i..(0.f)L«) s,., m,u..^ rAom 

G(p) = /P)((2.,0.0),0.«,(0,«)L..) - ^ = *("'(°- 

Neglecting again corrections of order l/L^, we have 

|J = 4^*K(0'^)^'*)' (A21) 

which proves Eq. (30). If we consider the corrections to scaling, this derivation shows that 
i?^ behaves essentially as U22 and C/4. The only difference is the presence of corrections 
due to the momentum-dependence of the scaling fields and to the specific definition of the 
correlation length: they scale as . . ., . . . Since u ^ 1^ m Eq. (30) they 

represent additional subleading corrections and can thus be neglected. This allows us to 
consider i?^, [/22, and on the same footing. 



APPENDIX B: SOME TECHNICAL DETAILS ON THE MC SIMULATIONS 

In our MC simulations we implement the standard Metropolis algorithm with a sequential 
update of the spins. We use a multispin^^ implementation, in which ribit = 64 systems are 
simulated in parallel. For each of them we use a different set of bonds {Jxy}- 

For the random numbers we use the SIMD-oriented Fast Mersenne Twister (SFMT)^''' 
generator. In particular, we use the genrand_res53() function that produces double-precision 
output. Independent random numbers are employed to generate the starting configurations 
for each disorder realization and in the parallel-tempering updates. In the latter case very 
few random numbers are used and thus, it takes virtually no extra time to use individual 
random numbers for each of the ^bit systems which are simulated in parallel. On the other 
hand, in order to save CPU time, we use the same sequence of random numbers for the 
local Metropolis update of any of the nbit systems. Though this choice does not lead to 
wrong estimates of the expectation values, it might create a statistical correlation among 
the ribit systems. However, since each of the ^bit systems correspond to a different set of 
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bond couplings Jxy we expect tliis effect to be negligible. Nevertheless, in order to ensure a 
correct estimate of the statistical error, in our jackknife analysis we put all ^bit systems that 
use the same sequence of random numbers in the same bin. In order to compute overlap 
observables, we performed runs for two systems with the same set {Jxy} in parallel. In our 
MC simulations a single Metropolis update of a single spin takes about 1.2 x 10"^ seconds 
on an Opteron CPU running at 2 GHz (this should be compared with the speed of the 
dedicated computer Janus,^^ the fastest computer simulating discrete spin models, which 
takes 2 x 10"^^ seconds to update an Ising spin). 

To reduce autocorrelations we used the random-exchange or parallel-tempering method. 
To this end, we divided the interval [/3min, /3max] into Np — 1 equal intervals A/?. The pa- 
rameter /9max was chosen such that ^(/3max)/-^ ~ 0.63 in most cases; in the latest runs we 
considered larger values, such that CiPmaxj/L ~ 0.66. The parameter /3min was chosen such 
that ^(/Sjnin) ^ L- We computed the observables in the neighborhood of /9max by using their 
second-order Taylor expansion around (3^ax- The coefficients of the expansion were estimated 
by measuring appropriate correlators with the energy. 

An elementary update unit consists in n^aet Metropolis sweeps followed by a replica ex- 
change of all pairs of systems at nearby temperatures. The different systems were sequentially 
visited, starting from those at /Smin and /?min + A/?. As a candidate for the exchange, we 
considered one of two replicas with equal probability. The acceptance probability for the 
exchange is min[l, exp(— A/3Aif)]. Since the measurement of the energy in our implemen- 
tation costs more CPU time than a Metropolis sweep, we chose ^met ^ 1 independent of 



The computation of disorder averages of products of thermal expectations requires par- 
ticular care. Indeed, naive estimators show a bias which may become larger than statistical 
errors. To avoid the problem we consider essentially bias-free estimators, defined following 
Rcf. 46. For this purpose wc divide the measurement phase of the run into 12 intervals. 
Between each pair of subsequent interval there is a dccorrelation phase. In total, the run 
consists of the following phases: Eq, Di Mi, D2 M2, -D12 M12. After some tests, we fixed 
the number of update steps for each of them. The equilibration phase Eq corresponds to 
20ntcmp elementary update units; the measurement phases correspond to ntcmp update 
units, while the length of Di is n^cxnp for i 7^ 7 and 5ntcmp for i = 7. Recall that each 
elementary update unit corresponds to nmet Metropolis sweeps of all systems and to one full 
tempering sweep. 

The presence of different measurement phases allows us to define bias-free quantities. To 
define [{A){B)] (for instance, this is relevant for the computation of U22) we average over 
the samples the quantity 



where fi{A)i is the average of estimates of A obtained in the measurement phase Mj. Analo- 
gously, to compute [{A){B){C)] (these correlators are necessary to compute the coefficients 
of the Taylor expansions around a given value of /3) , we average over the samples the quantity 
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TABLE IV: Summary of the parameters for the runs at p = 0.5. is the number of /3 vahics used 
in the parallel-tempering simulation. The rest of the notation is explained in the text. The CPU 
time refers to a single core of a dual core Opteron CPU running at 2.4 GHz. 
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0.58 


0.908 
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0.58 
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7 
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8 
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■J- \J \J \J \J \J 
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10 
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0.55 
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27 
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10 
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0.55 


0.896 


50 


11 


109779 


10 


300 


10 


0.54 
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12 
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10 
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10 


0.54 


8955 


308 


13 
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10 
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10 


0.54 


0.8955 
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14 
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10 


0.62 


0.8955 


361 


16 


24331 


10 
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20 


0.52 


0.895 


830 


20 


1542 


20 


2000 


32 


0.5125 


0.895 


658 


20 


2291 


50 


1500 


20 


0.625 


0.91 


1146 


24 


717 


25 
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32 


0.5125 


0.895 


826 


24 


1627 


50 


2000 


20 


0.625 


0.91 


1874 


28 


285 


60 


2500 


20 


0.6575 


0.895 
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In order to check equilibration, and decorrelation for the bias correction, we followed the 
suggestion of Ref. 30. We doubled the length of the run until the estimates of all observables 
were consistent within error bars. We performed this check only for the observables at /Saiax, 
because these are expected to be the most difficult ones for equilibration and decorrelation. 
Starting from disordered configurations, we determined the number of update steps nhaif 
that are needed to reach (averaged over samples) half of the equilibrium value of the overlap 
susceptibility. In total, the equilibration consisted of at least lOO^haif update steps. Using 
these methods to check equilibration, we came up with the choices summarized in Table IV. 
The parameters are not highly tuned, since we had the CPU time available on short notice. 
The runs that were done later have typically a larger /3min than those done earlier. The run 
for L = 28 is a bit at the edge of the criterion given above for equilibration. However, given 
the rather small number of samples (Ng — 18240), we are quite confident that the estimates 
are correct within the quoted error bars. 
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